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1 Introduction 

The classical description of the dynamics of a large set of neurons is based on de- 
terministic/stochastic differential systems for the excitatory-inhibitory neuron net- 
work [1, 2]. One of the most classical models is the so-called noisy leaky integrate 
and fire (NLIF) model. Here, the dynamical behavior of the ensemble of neurons is 
encoded in a stochastic differential equation for the evolution in time of membrane 
potential v(t ) of a typical neuron representative of the network. The neurons relax to- 
wards their resting potential Vl in the absence of any interaction. All the interactions 
of the neuron with the network are modelled by an incoming synaptic current I(t). 
More precisely, the evolution of the membrane potential follows, see [3-8] 



where C m is the capacitance of the membrane and gi is the leak conductance, nor- 
mally taken to be constants with r m = gt/C m »2ms being the typical relaxation 
time of the potential towards the leak reversal (resting) potential Vl ^ —70 mV. Here, 
the synaptic current takes the form of a stochastic process given by: 



where <5 is the Dirac Delta at 0. Here, Je and Jj are the strength of the synapses, Ce 
and Cj are the total number of presynaptic neurons and t' E j and t\ j are the times of the 
j th -spike coming from the ; th -presynaptic neuron for excitatory and inhibitory neu- 
rons respectively. The stochastic character is embedded in the distribution of the spike 
times of neurons. Actually, each neuron is assumed to spike according to a stationary 
Poisson process with constant probability of emitting a spike per unit time v. More- 
over, all these processes are assumed to be independent between neurons. With these 
assumptions the average value of the current and its variance are given by \±c — bv 
with b — CeJe — C\Ji and cr^ = (Cg/J + Cijf)v. We will say that the network is 
average-excitatory (average-inhibitory resp.) if b > 0 (b < 0 resp.). 

Being the discrete Poisson processes still very difficult to analyze, many authors in 
the literature [3-5, 7-9] have adopted the diffusion approximation where the synap- 
tic current is approximated by a continuous in time stochastic process of Ornstein- 
Uhlenbeck type with the same mean and variance as the Poissonian spike-train pro- 
cess. More precisely, we approximate I(t) in (1.2) as 



where B t is the standard Brownian motion, that is, B t are independent Gaussian pro- 
cesses of zero mean and unit standard deviation. We refer to the work [5] for a nice 
review and discussion of the diffusion approximation which becomes exact in the 
infinitely large network limit, if the synaptic efficacies Je and // are scaled appro- 
priately with the network sizes Ce and C\. 



C m — = -gL(V-V L ) + I(t), 



(1.1) 




(1.2) 



/ (t) dt ft* fji c dt + &c dB t , 
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Finally, another important ingredient in the modelling comes from the fact that 
neurons only fire when their voltage reaches certain threshold value called the thresh- 
old or firing voltage Vf w —50 mV. Once this voltage is attained, they discharge 
themselves, sending a spike signal over the network. We assume that they instanta- 
neously relax toward a reset value of the voltage Vr «s —60 mV. This is fundamental 
for the interactions with the network that may help increase their membrane potential 
up to the maximum level (excitatory synapses), or decrease it for inhibitory synapses. 
Choosing our voltage and time units in such a way that C m — gL — 1, we can sum- 
marize our approximation to the stochastic differential equation model (1.1) as the 
evolution given by 

dV = (-V + V L + ^c)dt + a c dB t (1.3) 

for V < Vp with the jump process: V(t£) = Vr whenever at to the voltage achieves 
the threshold value V(t~) = Vf', with Vl < Vr < Vf - Finally, we have to specify the 
probability of firing per unit time of the Poissonian spike train v. This is the so-called 
firing rate and it should be self-consistently computed from a fully coupled network 
together with some external stimuli. Therefore, the firing rate is computed as v — 
Vext + N(t), see [5] for instance, where N(t) is the mean firing rate of the network. 
The value of N{t) is then computed as the flux of neurons across the threshold or 
firing voltage Vf- We finally refer to [10] for a nice brief introduction to this subject. 

Coming back to the diffusion approximation in (1.3), we can write a partial dif- 
ferential equation for the evolution of the probability density p(v, t) > 0 of finding 
neurons at a voltage v e (— oo, Vf] at a time t > 0. A heuristic argument using Ito's 
rule [3-5, 7-9, 11] gives the backward Kolmogorov or Fokker-Planck equation with 
sources 

^(v,t) + ^-[h(v,N(t))p(v,t)]-a(N(t))^(v,t) 

dt dv dv z (1.4) 

= S(v-V R )N(t), v<V F , 

with h(v, N(t)) — —v + Vf + [i c and a(N) — er^ /2. We have the presence of a source 
term in the right-hand side due to all neurons that at time t > 0 fired, sent the signal 
on the network and then, their voltage was immediately reset to voltage Vr . More- 
over, no neuron should have the firing voltage due to the instantaneous discharge of 
the neurons to reset value Vr, then we complement (1.4) with Dirichlet and initial 
boundary conditions 

p(V F ,t)=0, p(-oo,0 = 0, p(v,O)=p°(v)>0. (1.5) 

Equation (1.4) should be the evolution of a probability density, therefore 

/Vf r Vf 

p(v,t)dv= / p (v)dv = l 
-oo J —oo 

for all t > 0. Formally, this conservation should come from integrating (1 .4) and using 
the boundary conditions (1.5). It is straightforward to check that this conservation for 
smooth solutions is equivalent to characterize the mean firing rate for the network 
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N(t) as the flux of neurons at the firing rate voltage. More precisely, the mean firing 
rate N(t) is implicitly given by 

N(t):=-a(N(t))^-(V F ,t)>0. (1.6) 

Here, the right-hand side is nonnegative since p > 0 over the interval [— oo, Vf] 
and thus, £(Vf, t) < 0. In particular this imposes a limitation on the growth of the 
function N \-> a(N) such that (1.6) has a unique solution N. Let us mention that a 
rigorous passage from the stochastic differential equation with jump processes (1.3) 
to the nonlinear equation (1 .4)-(l .6) is a very interesting issue but outside the scope of 
this paper, see related results in which nonlinearities are nonlocal functionals in [12, 
13]. 

The above Fokker-Planck equation has been widely used in neurosciences. Often 
the authors prefer to write it in an equivalent but less singular form. To avoid the 
Dirac delta in the right hand side, one can also set the same equation on (— oo, Vr) U 
(Vr, Vf] and introduce the jump condition 

9 9 
p(Yg, t) = p(V+, t), —piV^, t) - —p(V+, t) = N(t). 
dv dv 

This is completely transparent in our analysis which relates on a weak form that 
applies to both settings. 

Finally, let us choose a new voltage variable by translating it with the factor 
Vl + bvext while, for the sake of clarity, keeping the notation for the rest of values 
of the potentials involved Vr < Vf- In these new variables, the drift and diffusion 
coefficients are of the form 

h(v,N) = -v + bN, a(N) = a 0 + ai N, (1.7) 

where b > 0 for excitatory-average networks and b < 0 for inhibitory-average net- 
works, «o > 0 and a\ > 0. Some results in this work can be obtained for some more 
general drift and diffusion coefficients. The precise assumptions will be specified on 
each result. Periodic solutions have been numerically reported and analysed in the 
case of the Fokker-Planck equation for uncoupled neurons in [14, 15]. Also, they 
study the stationary solutions for fully coupled networks obtaining and solving nu- 
merically the implicit relation that the firing rate N has to satisfy, see Section 3 for 
more details. 

There are several other routes towards modeling of spiking neurons that are related 
to ours and that have been used in neurosciences, see [16]. Among them are the deter- 
ministic I&F models with adaptation which are known for fitting well experimental 
data [17]. General models of this type were unified and studied in terms of neuronal 
behaviors in [18]. In this case it is known that in the quadratic (or merely superlinear) 
case, the model can lead to blow-up [19] in the absence of a fixed threshold. We point 
out that the nature of this blow-up is completely different from the one discussed in 
this paper. One can also introduce gating variables in neuron networks and this leads 
to a kinetic equation, see [20] and the references therein. Another method consists 
in coding the information in the distribution of time elapsed between discharges [21, 
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22], this leads to nonlinear models that exhibit naturally periodic activity and blow-up 
cannot happen. Nonlinear IF models are able to produce different patterns of activity 
and excitability types while linear models do not. 

In this work we will analyse certain properties of the solutions to (1.4)-(1.5) with 
the nonlinear term due to the coupling of the mean firing rate given by (1.6). Next 
section is devoted to a finite time blow-up of weak solutions for (1.4)-(1.6). In short, 
we show that whenever the value of b > 0 is, we can find suitable initial data con- 
centrated enough at the firing rate such that the defined weak solutions do not exist 
for all times. We remark that, in the same sense that Brunei in [4], we use the term 
asynchronous for network states for which the firing rate tends asymptotically to con- 
stant in time, while we denote by synchronous those for which this does not happen. 
Therefore, a possible interpretation of the blow-up is that synchronization occurs in 
the model since the firing rate diverges for a fixed time creating possibly a strong par- 
tial synchronization, that is, a part of the network firing at the same time. Although 
one could also consider the blow-up as an artifact of these solutions, since neurons 
firing arbitrarily fast is not biologically plausible. As long as the solution exists in 
the sense specified in Section 2, we can get a priori estimates on the L^-norm of 
the firing rate. Section 3 deals with the stationary states of (1.4)-(1.6). We can show 
that there are unique stationary states for b < 0 and a constant but for b > 0 dif- 
ferent cases may happen: one, two or no stationary states depending on how large 
b is. In Section 4, we discuss the linear problem b = 0 with a constant for which 
the general relative entropy principle applies implying the exponential convergence 
towards equilibrium. Finally by means of numerical simulations, in Section 5 we il- 
lustrate the results of previous sections about blow-up and steady states. Moreover, 
this numerical analysis allows us to conjecture about nonlinear stability properties of 
the stationary states: in case of only one steady state it is asymptotically stable and 
in case of two different stationary solutions the results show that the one with lower 
firing rate is locally asymptotically stable while the one with higher stationary firing 
value is either unstable or with a very small region of attraction. Our results and sim- 
ulations describe situations which can be identified with neuronal phenomena such 
as synchronization/asynchronization of a network and bistability networks. Bi- and 
multi-stable networks have been used, for instance in models of visual perception 
and decision making [23-25]. Our analysis in Sections 2, 3 and 5 imply that this 
simple model encodes complicated dynamics, in the sense that, only in terms of the 
connectivity parameter b, very different situations can be described with this model: 
blow-up, no steady state, only one steady state and several stationary states. 



2 Finite time blow-up and a priori estimates for weak solutions 

Since we study a nonlinear version of the backward Kolmogorov or Fokker-Planck 
equation (1.4), we start with the notion of solution: 

Definition 2.1 We say that a pair of nonnegative functions (p, N) with p e 
L°°(R+; L l + (-oo, V F )), N e L] c+ (K + ) is a weak solution of (1.4)-(1.7) if 
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for any test function <j){v,t) e C°°((— oo, Vf] x [0, T]) such that 
L°°((-oo, Vf) x (0, T)), we have 



&± 30 
3t> 2 ' U 3« 



p{v,t) 



90 30,. ... d z <p 
9f 3n z 



dvdt 



= [ N(t)[<j>(V R ,t)-<j>(V F ,t)]dt (2.1) 
Jo 

p°(v)4>(0,v)dv- p(v,T)<j>(T,v)dv. 

-oo J—oo 



Here, the notation L p (£2), 1 < p < oo, refers to the space of functions such that 
f p is integrable in £2, while L°° corresponds to the space of bounded functions in £l. 
The set of infinitely differentiable functions in Q. is denoted by C°° (£2) used as test 
functions in the notion of weak solution. These non-negativity assumptions are rea- 
sonable. Indeed, for a given N(t), if we were to replace N by N + in the right hand 
side of (1.4), we obtain a linear equation which solution is non-negative and (1.6) 
gives N > 0; that this fixed point may work is a more involved issue, since we prove 
that there are not always global solutions, which requires functional spaces and this 
motives the a priori estimates that we derive at the end of this section. 

Let us remark that the growth condition on the test function together with the 
assumption (1.7) imply that the term involving h(v, N) makes sense. By choosing 
test functions of the form \[r(t)4>(v), this formulation is equivalent to say that for all 
(j>{v) e C°°((-oo, V F ]) such that i>§£ e L°°((-oo, V»), we have that 



d C Vf 
■ I (p(v)p(v, t)dv 

J—oo 



dt t 

(2.2) 

' -TACv.AO + a— j p(v,t)dv + N(t)[<t>(y R )-(t>(yF)] 

-oo \_OV OV J 



holds in the distributional sense. It is trivial to check that weak solutions conserve the 
mass of the initial data by choosing <f> — 1 in (2.2), and thus, 

( f p(v,t)dv= { F p°(v)dv=l. (2.3) 

J—oo J—oo 

The first result we show is that global-in-time weak solutions of (1.4)-(1.6) do not 
exist for all initial data in the case of an average-excitatory network. This result holds 
with less stringent hypotheses on the coefficients than in (1.7) with an analogous 
notion of weak solution as in Definition 2.1. 



Theorem 2.2 (Blow-up) Assume that the drift and diffusion coefficients satisfy 

h(v,N) + v>bN and a(N)>a m >0, (2.4) 

for all — oo < v < Vf and all N > 0, and let us consider the average-excitatory 
network where b > 0. Choose p. > max(— , V). If the initial data is concentrated 
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enough around v — Vf, in the sense that 

• V F 



/ 



e IJ - v p°(v)dv 



is close enough to e^ Vr , then there are no global-in-time weak solutions to (1.4)- 
(1.6). 

Proof We choose a multiplier 4>{v) — e^ v with pi > 0 and define the number 

, ct>(v F )-4>(v R ) 

X — > 0 

bp, 

by hypotheses. For a weak solution according to (2.1), we find from (2.2) that 
d C Vf 

— I (p(v)p(v,t)dv 

Ml J — OO 

rV F 

>p,l (bN(t)-v)<t>(v)p(v,t)dv 

J — OO 

(2.5) 

/v F 
(f>(v)p(v, t)dv - \bp,N(t) 
-OO 

fVp 

> /x / <p(v)p(v, t)dv [bN(t) + pta m - Vf] - Xp,bN(t), 

J —OO 

where (2.4) and the fact that v e (-co, Vp) was used. Let us now choose p large 
enough such that pa m — V F > 0 according to our hypotheses and denote 

• V F 



/vf 
<t>(v)p(v,t)dv, 
-OO 

;M„(t) > bnN(t)[M^(t) - A]. 



which satisfies 

d 

dt' 

If initially M M (0) > X and using Gronwall's Lemma since N(t) > 0, we have that 
M^(t) > X, for all t > 0, and back to (2.5) we find 



d f Vf f Vf 

— / <p(v)p(v, t)dv > p{p,a m — Vf ) / cf){v)p(v,t)dv 

dt J —OO J —OO 

which in turn implies, 

( * </>(v)p(v,t)dv>ei l(lJ - a >«- VF) ' ( " <j>(v)p°(v)dv. 

J — OO J —OO 

On the other hand, since p(v, t) is a probability density, see (2.3), and p > 0 then 

fVp 

/ (p(v)p(v,t)dv<e^ VF , 

J — CO 
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leading to a contradiction. 

It remains to show that the set of initial data satisfying the size condition (0) > 
X is not empty. To verify this, we can approximate as much as we want by smooth 
initial probability densities an initial Dirac mass at Vf which gives the condition 



e nV F >X — together with [ia m > Vf- 

bp, 



This can be equivalently written as 



1 _ e -ii(V F -V R ) y p 

p > and p > — . 

b a m 

Choosing a > max(f , — ), these conditions are obviously fulfilled. □ 

As usual for this type of blow-up result similar in spirit to the classical Keller- 
Segel model for chemotaxis [26, 27], the proof only ensures that solutions for those 
initial data do not exist beyond a finite maximal time of existence. It does not charac- 
terize the nature of the first singularity which occurs. It implies that either the decay at 
infinity is false, although not probable, implying that the time evolution of probability 
densities ceases to be tight, or the function N(t) may become a singular measure in 
finite time instead of being an Lj oc (R + ) function. Actually, in the numerical com- 
putations shown in Section 5, we observe a blow-up in the value of the mean firing 
rate in finite time. To continue the solution will need a modification of the notion of 
solution introduced in Definition 2.1. It would be useful, since the firing rate does 
not become constant in time and consequently a possible interpretation is that syn- 
chronization occurs, a phenomena that is of interest to neuroscientists, see also the 
comments in the introduction and the conclusions. 

Although in this paper the nature of blow-up is not mathematically identified, we 
devote the rest of the section to prove some a priori estimates which shed some light 
on this direction. To be more precise, our estimates indicate that this blow-up should 
not come from a loss of mass at v & — oo, or a lack of fast decay rate because the 
second moment in v is controlled uniformly in blow-up situations. We obtain these 
a priori bounds with the help of appropriate choices of the test function (p in (2.1). 
Some of these choices are not allowed due to the growth at — oo of the test functions. 
We will say that a weak solution is fast-decaying at — oo if they are weak solutions 
in the sense of Definition 2.1 and the weak formulation in (2.2) holds for all test 
functions growing algebraically in v. 

Lemma 2.3 (A priori estimates) Assume (1.7) on the drift and diffusion coefficients 
and that (p,N) is a global-in-time solution of (1.4)-(1.6) in the sense of Defini- 
tion 2.1 fast decaying at — oo, then the following a priori estimates hold for all 
T > 0: 

(i) Ifb> V F - V R ,then 

/V F I rVp \ 

{V F - v)p(v, t)dv < max I V F , / (V F - v)p°(v)dv J, 
-oo \ J —oo / 
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) f N(t)dt < V F T + f 
Jo J-t 



(b-V F +V R ) I N(t)dt <V F T + I (V F -v)p°(v)dv. 



(ii) Ifb < V F - Vr then 

• Vf / rV F 



/vf / rv F \ 

(V F -v)p(v,t)dv>minl V F , / (Vf — v)p (v)dv J. 
-00 \ J — oo J 



Moreover, if in addition a is constant then 
/o 



f 

Jo 



N{t)dt <(l + T)C(V F ,V R ,a, p°). (2.6) 



Proof Using (1.7) together with our decay assumption at — oo, we may use the test 
function (p (v) — V F — v > 0. Then (2.2) gives 



d [ Vf 

I cp(v)p(v,t)dv 

J —c 



dt j 



[v - bN(t)]p(v, t) dv + N(t)(V F - V R ). 

-oo 



This is also written as 

— / (p(v)p(v, t)dv + I 4>(v)p(v,t)dv 

at J-oo J-oo (2.7) 

= V F -N(t)[b-(V F -V R )]. 



To prove (i), we notice that with our condition on b, the term in N(t) is nonpositive 
and the first result follows from Gronwall's inequality. The second result just follows 
after integration in time. 

To prove (ii), we first use again (2.7) and, because the term in N(t) is non- 
negative, we find the first result. To obtain the second estimate (2.6), given e e 
(0, (V F — Vr)/2), we can always choose a smooth truncation function 4> € {v) € C 2 
such that 



<t> e (V F ) = l, Mv) = 0 foTv<V R , (p' € (v)>0, 4>"(v)>0, 
with <j)"(v) — 0 outside the interval (Vr, Vr + e) such that 

»>)->- ., T7 focaHve(Vn,V F ], as e -► 0, (2.8) 
V F - Vr 

and thus (/>" e L°°(— oo, V F ) with size of order of e~ l . In other words, we have 
chosen a C 2 uniform approximation of the truncation (v — Vr)+/{V f — Vr) with 
x+ — max(x,0) obtained by integrating twice a smooth suitable approximation of 
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the S(V - V R )/(V F - Vr). Then, equation (2.2) gives 
d C Vf 

— I <p € (v)p(v,t)dv + N(t) 

at j Vr 

= I (p' € (v)(-v + bN(t))p(v,t)dv + a 4>"(v)p(v,t)dv. 
JVr Jv r 

Since 4>' € {v) is positive and non-decreasing and using (2.8), we get for all 0 < e < 1 
there exists eo small enough, such that for all 0 < e < €q: 



/ <j>' € (v)(-v + bN{t))p{v,t)dv 
Jv R 



' v F 

<(/>' e (VF)(max(\V F \,\V R \)+bN(t)) (2.9) 

< (l + e)(max(\V F \,\V R \) + bN(t)) 
V F -V R 

Due to the hypotheses b < V F — Vr, for any y such that 0 < y < 1 — b/(V F — V R ), 
we can find e small enough such that 

(l + ~s)b 

0 < v < 1 . 

V F - Vr 

Taking into account (2.8), (2.9), and that p(v, t) is a probability density, we have for 
0 < € < €q small enough 



— / 4> € {v)p(y,t)dv + yN(t)< — + a\\^\\ L ^v R ,v P y 

dt Jy R V F - Vr 



Choosing now e = eo/2 for instance, integration in time of the last inequality leads 
to the desired inequality (2.6). □ 

Corollary 2.4 Under the assumptions of Lemma 2.3 and assuming v 2 p°(v) e 
h} (— oo, Vp ) and 0 < b < Vp — Vr, then the following a priori estimates hold: 

(i) If additionally a is constant, for all t > 0 we have 

•V F 



/ v 2 p(v,t)dv<C(l+t). 

J—oo 



(ii) If additionally -bmin(V F , f^(V F - v)p°(v)dv) + a\ + bVp + %^ ^ °> 
then 

/Vp / /■ v e 

v 2 p(v, t)dv < max I ao, I v 2 p°(v,t)dv 
-oo \ J—oo 
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Proof We use again the weak formulation (2.1) with 4>{v) — v /2 as test function 
and get 



dt 



[ VF V Z f F 7 

/ — p(v,t)dv+ / v p{v,t)dv 

J —oo 2 J —oo 

/ vp(y,f)dv + a(N(t)) + 

J —oo 



= bN(t) 



= bN(t) 



N(ty 



vi - vi 



/ (v 

J—oo 



<a 0 + N(t) 



- V F )p(v,t)dv + a(N(tj) + N(t) 



bV F 



Vi-V 



2 -\ 



b min I Vf 



f 

J— ( 



(V F -v)p u (v)dv )+a l +bV F 



Vi-V, 



thanks to the first statement of Lemma 2.3(ii). 

To prove (i), we just use the second statement of Lemma 2.3(ii) valid for a constant 
which tells us that the time integration of the right-hand side grows at most linearly 
in time and so does f^!L v 2 p(v, t)dv. 

To prove (ii), we just use that the bracket is nonpositive and the results follows. □ 



3 Steady states 

3.1 Generalities 

This section is devoted to find all smooth stationary solutions of the problem (1.4)- 
(1.6) in the particular relevant case of a drift of the form h(y) — Vq(N) — v. Let us 
search for continuous stationary solutions p of (1.4) such that p is C 1 regular except 
possibly at V — Vr where it is Lipschitz. Using the definition in (2.2), we are then 
allowed by a direct integration by parts in the second derivative term of p to deduce 
that p satisfies 



9 

9 v 



(v - V 0 (N))p(v) + a(N)^-p(v) + NH(v - V R ) 



= 0 (3.1) 



in the sense of distributions, with H being the Heaviside function, that is, H(u) — 1 
for u > 0 and H{u) — 0 for u < 0. Therefore, we conclude that 

(v - V Q (N))p + a(N)^- + NH(v - V R ) = C. 
dv 

The definition of N in (1.6) and the Dirichlet boundary condition (1.5) imply C — 0 
by evaluating this expression at v — Vf- Using again the boundary condition (1.5), 
P(Vf) = 0, we may finally integrate again and find that 

N (v-V 0 (N)) 2 fV F (w-V 0 (N)) 2 

p{v) = e / e 2 "W> H[w-V R ]dw 

a(N) J v 
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which can be rewritten, using the expression of the Heaviside function, as 



p(v) ■■ 



N (v-vq(n)) 2 r v F 



a(N) 



L 



e wm dw. 



(3.2) 



max(«, Vr) 



Moreover, the firing rate in the stationary state N is determined by the normalization 
condition (2.3), or equivalently, 



N 



V F 



{v-Vq(N)) 2 fVp (w-V 0 (N)) z 

e WW / e WW) dw 

J max(u, Vr) 



dv. 



(3.3) 



Summarizing, all solutions p of the stationary problem (3.1), with the above referred 
regularity, are of the form given by the expression (3.2), where N is any positive 
solution of the implicit equation (3.3). 

Let us first comment that in the linear case Vo(N) = 0 and a(N) — «o > 0, we 
then get a unique stationary state given by the Dawson function 



Poo(v) = e 



V 2 fV F 

la o J e 2a o di 

J max(v,Vf() 



(3.4) 



with jVoo the normalizing constant to unit mass over the interval (— oo, Vp]. 

The rest of this section is devoted to find conditions on the parameters of the model 
clarifying the number of solutions to (3.3). With this aim, it is convenient to perform 
a change of variables, and use new notations 



z — 



v-Vq 
-Ja 



w-Vo 



V F -Vq Vr- V 0 
Wf— }= , WR — }= ■ 



(3.5) 



where the dependency has been avoided to simplify notation. Then, as in [3], we 
can rewrite the previous integral (and thus the condition for a steady state) as 



1 

N ~ 
I(N) 



I(N), 



/ w f r z 2 ru>F 
-oo L ^max(z 



e 2 du 



max(z, wr) 



(3.6) 



dz. 



Another alternative form of I(N) follows from the change of variables s — (z — u)/2 
and s — (z + u)/2 to get 



I(N)=2 



0 rWf-\-s 



/ 7 



e~ 2s 'dSds. 



CO J WR-\-S 



f 

J — oc 



-2s Wf 



— e 



-2swr 



)ds, 



and consequently, 



I(N) 



I 



oo e -s A /2 



SWf sw 



R )ds. 



(3.7) 
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3.2 Case of a (AO = ao 

We are now ready to state our main result on steady states. 

Theorem 3.1 Assume h(v, N) — bN — v, a(N) — ao is constant and Vo = bN . 

(i) For b < 0 and b > 0 small enough there is a unique steady state to (1.4)-(1.6). 

(ii) Under either the condition 

0<b <V F -V R , (3.8) 

or the condition 

Q<2a 0 b<(V F -V R ) 2 VR, (3.9) 

then there exists at least one steady state solution to (1.4)-(1.6). 

(iii) If both (3.9) and b > V F — Vr hold, then there are at least two steady states to 
(1.4>(1.6). 

(iv) There is no steady state to (1.4)-(1.6) under the high connectivity condition 

Z?>max(2(V F - V R ),2V F I (0)). (3.10) 



Remark 3.2 It is natural to relate the absence of steady state for b large with blow- 
up of solutions. However, Theorem 2.2 in Section 2 shows this is not the only possible 
cause since the blow-up can happen for initial data concentrated enough around Vf 
independently of the value ofb>0. See also Section 5 for related numerical results. 

Proof Let us first study properties of the function I(N). To do that, we rewrite (3.7) 
as 

sVp sV R 



I(N)= e~ s /2 e ds. 

Jo s 



sVp sV R 

Taking the function f{s) — e^> — e^o and Taylor expanding up to second order 
at s = 0, we get /<» - /(0) - f'(0)s = f"(6)s 2 /2 with /(0) = 0, /'(0) = (V F - 
Vr)/^/oo, and 6 e (0, s). It is easy to see that 

/ v I b _Xf_ Vl °-Xk 
\f"(9)\ <max -J-e^o ,-^-e^o 
\2ao 2ciq 

for all 0 e (0, s). By distinguishing the cases based on the signs of Vp and Vr, this 
Taylor expansion implies that 



V F -V R 



max(V|, Vl) ,aa ^ '' A' 11 



2 «o (3.11) 



:= C 0 se v^o 
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for all s > 0. Then, a direct application of the dominated convergence theorem and 
continuity theorems of integrals with respect to parameters show that the function 
I(N) is continuous on N on [0, oo). Moreover, the function I(N) is C°° on N since 
all their derivatives can be computed by differentiating under the integral sign by 
direct application of dominated convergence theorems and differentiation theorems 
of integrals with respect to parameters. In particular, 



/'(AO = 
and for all integers k > 1 , 

/(*)(#) = (-!)* 



b 
ao 



b 



f 

Jo 



s 2 /2l sw F 



e r — e 



SW R 



)ds, 



f 

Jo 



<)ds. 



As a consequence, we deduce: 

1. Case b < 0: I{N) is an increasing strictly convex function and thus 

lim I{N) — oo. 

2. Case b > 0: I(N) is a decreasing convex function. Also, it is obvious from the 
previous expansion (3.11) and dominated convergence theorem that 

lim I(N) = 0. 



It is also useful to keep in mind that, thanks to the form of I(N) in (3.6), 

7(0) < y/2n[w F (0) - w R (0)]e max(w > ) ' w2 F m/2 

max(y2 ; y2) 



'2tt — — exp 



Now, let us show that for b > 0, we have 



lim NI(N) : 

N—>oo 



2a 0 



Vf-Vr 



(3.12) 



< oo. 



(3.13) 



Using (3.11), we deduce 



NI(N)-N 



Vf-Vr 



«o Jo 



f 

Jo 



s I 2 



e ^ ds 



< CqN / se 



I 

Jo 



sbN Jmax(|V F |,|V^|) 



Me 



ds. 



A direct application of dominated convergence theorem shows that the right hand 



side converges to 0 as N —>■ oo since s A 7 exp(— is a bounded function uniform 
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in N and s. Thus, the computation of the limit is reduced to show 

f 

Jo 



lim N I e - s2/2 ~^ds= ^ 
N-*-oo 



(3.14) 



With this aim, we rewrite the integral in terms of the complementary error function 
defined as 



erfc(x) :- 



) poo 

= e ~'~ 

* J.x 



and then 



L 



00 - s 2 /2 -ibM- 

e - M 



f 

Jo 



■ Jn ^ ( bN \ 

ds — -^—e^o erfc ; . 

V2 WW 



Finally, we can obtain the limit (3.14) using L'Hopital's rule 



11 

JO 



lim N I e ' ' M 



r- erfc(-^=) 
J*ds=^= lim ^° 



V2 N- 



b 2 N 2 
2 "0 



2«2 



~Jl lim 



/2«o 



«0 



b 2 



b 2 N 2 



b 2 N 2 

-e ^0 



A* 2 



With this analysis of the function / (AO we can now proof each of the statements of 
Theorem 3.1: 



Proof of (i) Let us start with the case b < 0. Here, the function /(AO is increas- 
ing, starting at 7(0) < 00 due to (3.12) and such that 

lim /(AO = 00. 

A?->oo 

Therefore, it crosses to the function 1/N at a single point. 

Now, for the case b > 0 small, we first remark that similar dominated convergence 
arguments as above show that both / (AO and /'(AO are smooth functions of b. More- 
over, it is simple to realize that / (AO is a decreasing function of the parameter b. Now, 
choosing 0 < b < b* < ( V F - V R )/2, then I(N) > I*(N) for all N > 0 where I*(N) 
denotes the function associated to the parameter b*. Using the limit (3.13), we can 
now infer the existence of A* > 0 depending only on b* such that 

Vf - Vr 

A?/ (AO > Nh{N) > " > 1 

2b* 

for all N > AO. Therefore, by continuity of NI(N) there are solutions to NI(N) = 1 
and all possible crossings of I(N) and l/N are on the interval [0, N*]. We observe 
that both /(AO and /'(AO converge towards the constant function 1(0) > 0 and to 
0 respectively, uniformly in the interval [0, A/*] as b —> 0. Therefore, for b small 
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A/7 (AO is strictly increasing on the interval [0, A 7 *] and there is a unique solution to 
NI(N) = l. a 

Proof of (ii) 

Case of (3.8) The claim that there are solutions to NI (N) — 1 for 0 < b < V F - Vr 
is a direct consequence of the continuity of I(N), (3.12) and (3.13). 

Case of (3.9) We are going to prove that I(N) > l/N for (v ^° Vf)2 < N < , 
which concludes the existence of a steady state since 7(0) < oo due to (3.12) implies 
that I(N) < l/N for small N. Condition (3.9) only asserts that this interval for N is 
not empty. To do so, we show that 



r^n ( V R- V F) 2 f „ 
I (N) > for N e 

2a 0 



Vr 



which obviously concludes the desired inequality / (AO > l/N for the interval of N 
under consideration. 

The condition ^ > N is equivalent to wr > 0, therefore, using (3.5) and the 
expression for I(N) in (3.6), we deduce 

"Wf 



I(N)> 



2 fWp 



J ^ 



du 



dz > 



f WF [ z 2 f 



2 ru>F u 2 
e~*T I du 



dz- 



Since z > 0 and e 2 is an increasing function for u > 0, then e 2 > e 2 on [z, u>f], 
and we conclude 



/(A0> 



I I du 



dz ■ 



(Vr - VfY 
2fl 0 



□ 



Proof of (Hi) Under the condition (3.9), we have shown in the previous point the 
existence of an interval where I(N) > I/N. On one hand, 1(0) < 00 in (3.12) im- 
plies that I(N) < l/N for N small and the condition b > Vf — Vr implies that 
I(N) < I/N for N large enough due to the limit (3.13), thus there are at least two 
crossings between I(N) and l/N. □ 

Proof of (iv) Under assumption (3.10) for b, it is easy to check that the following 
inequalities hold 

1(0) < l/N foxN<2V F /b (3.15) 

and 

y F _ v » 1 

— ^<— fovN>2V F /b. (3.16) 

bN — Vf N ' 
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We consider N such that /V > Vf /b, this means that wp < 0. We use the formula 
(3.7) for I(N) and write the inequalities 



oo „ poo 

e -(s-w F ) 2 /2 



I(N) < (w F - wr) [ e- s2/2 e SWF = (w F - w R )e w F^ 2 f 
Jo Jo 

= (w F - w R )e w2 F> 2 r e~ s2 > 2 < (w F - w R )e w2 F' 2 f 

J—WF J- 



oo 

_c2 / 



s z /2 

v F -v R 



Va 0 \w F \ 

where the mean- value theorem and w F < 0 were used. Then, we conclude that 

rr*n V F~V R V F - V R 

I(N) < — = — — — , for N > V F b. 

*Ja Q \w F | bN — V F 

Therefore, using Inequality (3.16): 

1 

I(N)<—, foTN>2V F /b 

and due to the fact that / is decreasing and Inequality (3.15), we have I(N) < 1(0) < 
l/N, for /V < 2Vp/b. In this way, we have shown that for all N, I(N) < l/N and 
consequently there is no steady state. □ 



Remark 3.3 The functions I (N) and 1 /N are depicted in Figure 1 for the case 
Vq{N) — bN and a(N) — ao illustrating the main result: steady states exist for small 
b and do not exist for large b while there is an intermediate range of existence of 
two stationary states. The numerical plots of the function N I (N) might indicate that 
there are only three possibilities: one stationary state, two stationary states and no 
stationary state. However, we are not able to prove or disprove the uniqueness of a 
maximum for the function N I (N) eventually giving this sharp result. 
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Remark 3.4 The condition (3.9) is not optimal and it can be improved by using 
one more term in the series expansion of the exponentials inside the integral of the 
expression of I(N) in (3.7). More precisely, if w F > wr > 0, we use 



e SWF ~ ^ = E^7 K "0>E^7 W " O ■ 
In this way, we get 



n=0 



n=0 



Jo \ v fl o 2 

_ (V f -V r )(V2^+(V f -Vr)) 



2a 0 



since 



e~ s / 2 ds = J-, I e~ s/2 sds = l and V 0 < Vr. 



f> CO „ / _, nC 

Jo V 2 Jo 

77ien, condition (3.9) caw fee improved to 

2a 0 b < V R (V F - V R )(^2^+ (V f - Vr)). 
Of course, this last inequality is not optimal either for the same reason as before. 

3.3 Case of a{N) — ciq + a\N 

We now treat the case of a(N) — «o + a\N, with ao, a\ > 0 with b > 0. Proceeding 
as above we can obtain from (3.7) the expression of its derivative 



/'(AO = 



+ 



dN 
d 



Vo(N) 
1 



dN\^ajN) 



(h(N) - I 2 (N)) 
{V F h(N)-V R I 2 (N)), 



(3.17) 



where 



f°° 1 f°° 7 

h(N)= e~ s/2 e SWF ds and h(N)= e~ s /2 e SWR ds. 
Jo Jo 



Therefore / (AO is decreasing since 



d 

dN 



Vb(A0 



2bao + ba\ N 



( 



dN\^/ajN)J (a 0 + ai N)y 2 



> 0 and 



<0. 
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Moreover, we can check that the computation of the limit (3.13) still holds. Actually, 
we have 



lim 



NI(N) 
V F -V R 



lim 



N 
Ja 



f 

Jo 



"■' 2 /2, 



erfc(-^L) 

e 2<* 
N 



erfc(to) 

= hm V 77 n~ — nm V 77 " 



— ^be 
Jit 



-b 2 a 



-2b 2 a 2 e- b2c,2 -e- b2a2 



where we have used the change a — -^|= and L'Hopital's rule. In the case b < 0, we 
can observe again by the same proof as before that I{N) — > oo when N — ► oo, and 
thus, by continuity there is at least one solution to NI(N) = 1. Nevertheless, it seems 
difficult to clarify perfectly the number of solutions due to the competing monotone 
functions in (3.17). 

The generalization of part of Theorem 3.1 is contained in the following result. We 
will skip its proof since it essentially follows the same steps as before with the new 
ingredients just mentioned. 

Corollary 3.5 Assume h(v,N) = bN — v, a(N) — ao + a\N with oq, a\ > 0. 

(i) Under either the condition b < Vf — Vr, or the conditions b > 0 and 2a^b + 
2a\ Vr < (Vf — Vr) 2 Vr, then there exists at least one steady state solution to 
(1.4>(1.6). 

(ii) If both 2aob + 2a\VR < (Vp — Vr) 2 Vr and b > Vp — Vr hold, then there are 
at least two steady states to (1.4)-(1.6). 

(iii) There is no steady state to (1.4)-(1.6) for b > max(2(VV — Vr), 2VfI(0)). 

These behaviours are depicted in Figure 2. Let us point out that if a is linear and 
b < 0, I (N) has not to be strictly increasing as in the constant diffusion case and it 
may have a minimum for N > 0. 




Fig. 2 Left figure: The function NI(N) against the constant 1 when a(N) is linear. For /:> = 0.5 we have 
considered a (N) = 0.5 + N/8, for b = 1.2: a(N) = 0.4+ AT/100 and for b = 8: a(N) = 6 + ^/100. Right 
figure: The function N I(N) against the constant 1 with b < 0 and a(N) = 1 + N. Here Vr = 1; Vf = 2. 
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At this point, a natural question is what happens with the stability of these steady 
states. In the next section we study it in the linear case, when the model presents 
only one steady state. An extension of the same techniques, entropy methods, to the 
nonlinear case is not straightforward at all. However, the results obtained in the linear 
case let us expect that for small connectivity parameter b the only steady state could 
be stable. On the other hand, numerical results presented in Section 5 give some 
numerical evidence of the stability/instability in different situations described by this 
model: only one steady state or two steady states, see that section for details. 



4 Linear equation and relaxation 

We study specifically the case of a linear equation that is b — 0 and a(N) = ciq, that 
is. 



dp(v,t) 9 



dt 



- ^-[vp(v, t)] - a 0 —^p(v, t) =S(v - V R )N(t), v < V F , 
dv A "- L 



dv 



p(V F ,t) = 0, N(t):=-ao—p(V F ,t)>0, a 0 > 0, 

OV 

p(v,0) = p°(v)>0, / p°(v)dv=l. 

J—oo 



(4.1) 



For later purposes, we remind that the steady state poo(v) given in (3.4) satisfies, for 
the case at hand, the equation 

d 9 2 
~^[vpoo(v)] - ao—^pooiv) =S(v - V^Noo, v < V F , 
dv dv A 



Poo(V F ) = 0, 

/ Poo(.v)dv=\ 
J—oo 



Noc--=-a 0 — Poo(V F )>0, 
dv 



(4.2) 



We will assume in this section that the initial data satisfies, for some C° > 0 



p°(v)<C 0 Poo (v). 



(4.3) 



Then, we take for granted that solutions of the linear problem exist, with the regularity 
needed in each result below and such that for all t > 0, 



p(v,t)<C {, p 00 (v). 



(4.4) 



The estimate (4.4) follows a posteriori from the relative entropy inequality that we 
state below, see more comments at the end of this section. This is an indication that the 
hypothesis for the initial data (4.3) will easily be propagated in time giving (4.4) with 
a well-posedness theory of classical fast-decaying solutions at hand. These solutions 
to (4.1) and (4.2) might be obtained by the method developed in [28] and will be 
analysed elsewhere. 
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We prove that the solutions to (4.1) converge in large times to the unique steady 
state poo(v). 

Theorem 4.1 (Exponential decay) Fast-decaying solutions to the equation (4.1) ver- 
ifying (4.4) satisfy 

/ Poc(v)\ — dv<e zfl ° vr / Poc(v)\ — dv. 

J-oo \ Poo(V) J J-oo \ Poo(v) J 



This result shows that no synchronization of neuronal activity can be expected 
when the network is not connected, since solutions tend to produce a constant firing 
rate, a very intuitive conclusion. Because the rate of decay is exponential, we also 
expect that small connectivity cannot create synchronization either, again an intuitive 
conclusion proved rigorously for the elapsed time structured model in [22]. Also the 
proof shows that two relaxation processes are involved in this effect: dissipation by 
the diffusion term and dissipation by the firing term. These relaxation effects are 
stated in the following theorem which also gives the natural bounds for the solutions 
to equation (4.1) (choosing G(u) — u 2 gives the natural energy space of the system, 
a weighted L 2 space). 

Theorem 4.2 (Relative entropy inequality) Fast-decaying solutions to equation (4.1) 
verifying (4.4) satisfy, for any smooth convex function G : K + — > K, the inequality 



d f v " fp(v,t)\ 

/ Poo(u)G(^^ )dv = -D G [p](t)<0, 

J — oo 



dl 



Poo(v) 



(4.5) 



with 

D G [p](t)=N 00 



Noo ) \Poo(v) 



N(t) p(v,t)\ G ,( p(v,t) 



Poo(v) 



Poo(v)J_ 



Vr 



V F 



a 0 I Pco(v)G 

-oo \Poo(v) 



p(v,t) 



3 (p(y,t) 



.dv \pco(v) 



dv>0. 



Proof of Theorem 4.1 The proof is standard in Fokker-Planck theory and follows by 
applying the relative entropy inequality (4.5) with G(x) — (x — l) 2 . Then, we obtain 

d [ Vf 7 pM A 2 , ^ , [ v * [ 3 fp(vj) \] 2 
-r Poo(v)\ — - - 1 ) dv < -2ao / Poo(v) —\ — — 1 I dv. 

dt J-oo \Poo(V) ) J-oo L3u\PooW / 

(4.6) 

To proceed further, we need an additional technical ingredient that we state and whose 
proof is postponed to an Appendix. 
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Proposition 4.3 There exists v > 0 such that 



[ Vf , J ?(«) 

J — OO 



Poo(v) 



pV F 

iv< Poo(v) 
J—oo 



3 / q(v) 



dv \poo(v) 



dv 



for all functions q such — e H l (poo(v) dv) and f_£g q(v) dv — 0. 



Poo 



Poincare's-like inequality in Proposition 4.3 applied to q — p — poo bounds the 
right hand side on (4.6) 



d [ v ? ,f p(v,t) 

I Poo(v)\ — - 1 ) dv<-2a 0 



dt J-oo ' ' ' \poo(v) 

Finally, the Gronwall lemma directly gives the result. 



V / Poo(v)\ — - 1 I dv. 

J—oo 



Poc(v) 



□ 



To show Theorem 4.2, which was used in the proof of Theorem 4.1, we need the 
following preliminary computations. 



Lemma 4.4 Given p a fast-decaying solution of (4.1) verifying (4.4), poo given by 
(3.4) and G a C 2 convex function, then the following relations hold: 



1JL 

dt poo 



2ao 3 \ d p d 2 p 

v-\ — Poo h a 0 —^ 

Poo dv J dv poo OV z poo 



S(v - V R )[ — 

Poo \Noo Poo 



(4.7) 



dt \poo 



2a 0 3 \ 3 / p 

Poo I — O I 

Poo dv J dv \poo 



3 2 r ( P 

dv L 



-a 0 G" 



Poo 



dv poo 



Noo ( N 

— S(v-V R ) I — 

Poo \Noo 



Poo 

Poo J \ Poo 



(4.8) 



9 J P 
w:PooG[ — 

at \Poo 



d 

dv 



VpooG 



Gil 

-NooHv-V R ) 



Poo 
2 



ao 



dv 2 



PooG 



Poo/\dv poo 

N p 



Proof Since f (-E-) — — 5^ 'i- we 

■> du v f>oo poo dv dv 



G\ 

Noo Poo / \ Poo 
l_ 3p_ J^SPoo „„, Qbtain 



Poo 



P dp c 



dp _ _9_/_P_ 
3d ^dvypooj Poo dv 



(4.9) 
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and 



d 2 p 
dv 2 



; Poo: 



d ( P \ dPoo , P 3 Poo 



dv 2 \pooj dv\pooj dv poo dv 2 



Using these two expressions in 



dt \poc 



1 dp 

Poo dt 
1 

Poo 



d r n d 2 1 

<5(u = V«)JV(0 + — [vp(v, f)J + a 0 —jp(v, t) 
dv dv 1 j 



we obtain (4.7). 

Equation (4.8) is a consequence of Equation (4.7) and the following expressions 
for the partial derivatives of G(— ): 

Poo 



dt \poo 



■ G ' [ -)U-)> M- 

^Poo/dtypooJ dv \poo 



g' 



Poo / dv \poo 



and 



dv 2 \poo 



Poo/\dv\p 



G' 



Poo 



dv 2 



Poo 



Finally, Equation (4.9) is obtained using Equation (4.8) and the fact that p^ is solu- 
tion of (4.2). □ 

Proof of Theorem 4.2 We integrate from — oo to Vp — a in (4.9) and let a tend to 0 + 
and use L'Hopital's rule 



r P(V,t) .. 8 £(v,t) N(t) 

hm = hm = . 

v^V F Poo (v) v^V F °P°°( V ) Noo 

Since p(v,t) < C°poo(v), then 
d rV F 



(4.10) 



dl 



/ PooG\ )dv-a 0 

J — oo 



Poo 



8 v 



PooG 



Poo 



V F 



-a 0 j_ 



PooG"[ -^-|| -I dv 



3 p 



■[(: 



N 



Poo/\dv poo 

p \g'( '" 



.Noo Poo J \Poo 

The Dirichlet boundary condition (1.5) implies that 



Poo 



Vr 



—ao 



PooG 



Poo 



V F 



dPoo n ( P 
— «o G\ 

dv \poo 



v F \NooJ' 



^ Springer 



Page 24 of 33 



Caceres et al. 



where we used that 



3 r ( P 

OV \poo 



= PooG = 

V F \PocJ PioOQ 



Poo/ \ ao ao Poo 



V F 

= 0, 

V F 



due to (4.10). Collecting all terms leads to the desired inequality. □ 

Let us finally remark that as a usual consequence of the General Relative Entropy 
principle (GRE) [29], the estimate (4.4) follows by choosing the convex function 
G(x) — (x — C°) A + . This shows that the bound (4.4) can be proved using (4.3) together 
with a well-posedness theory of classical fast-decaying at — oo solutions to (4.1). 



5 Numerical results 

We consider two different explicit methods to simulate the NNLIF (1 .4). The first one 
is based on standard shock-capturing methods for the advection term and standard 
centered finite differences for the second-order term. More precisely, the first order 
term is approximated by finite difference WENO-schemes [30]. 

The second numerical method is based on another finite difference scheme for 
the Fokker-Planck equation proposed in the literature called the Chang-Cooper 
method [31]. This method was also used in [20] for a computational neuroscience 
model with variable voltage and conductance. In order to use this method, the first 
step is to rewrite the Fokker-Planck equation (1.4) in terms of the Maxwellian 

(v~bN) 2 

M(v) — e W as follows, 



^(v,t)-a(N(t))^- 
dt OV 



3 / p(v, t) 
M(v,t) ' ' 



dv \M(v,t) 



= N(t)S(v-V R ). 



Then, the Chang-Cooper method performs a kind of 8 -finite difference approximation 
of p/M, see [31] for details. The Chang-Cooper method presents difficulties when 
the firing rate becomes large and the diffusion coefficient a (N) is constant. More pre- 
cisely, given a(N) = ao and b > 0, if TV is large, the drift of the Maxwellian, in terms 
of which is rewritten the Fokker-Planck equation, practically vanishes on the interval 
(— oo, Vf] and this particular Chang-Cooper method is not suitable. Whenever a(N) 
is not constant, this problem disappears. 

Summarizing, we consider two different schemes for our simulations: the first 
one is based on WENO-finite differences as described above, and the second one 
by means of the cited Chang-Cooper method. In both cases the evolution on time is 
performed with a TVD Runge-Kutta scheme. In Section 2 of [20] these schemes are 
explained in details and we refer to [30, 31] for a deeper analysis of them. 

In our simulations we consider a uniform mesh in v, for v e [V m ,„, Vf]. The 
value Vmin (less than Vr) is adjusted in the numerical experiments to fulfill that 
p(V,nin, t) & 0, while Vf is fixed to 2 and Vr = 1. As initial data we have taken 
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two different types of functions: 

• Maxwellians: 

p°(v)= e 2 "o , (5.1) 

where the mean Do and the variance are chosen according to the analyzed phe- 
nomenon. 

• Stationary Profiles (3.2) given by 

N (v-V 0 (N)) 2 fV F (w-V 0 (N)) 2 

P(V) - e WN) / e 2a(N) dw, 

a(N) Jmax(v,V R ) 

with N an approximate value of the stationary firing rate. We typically consider 
this kind of initial data to analyze local stability of steady states. 

Steady states. - As we show in Section 3, for b positive there is a range of values 
for which there are either one or two or no steady states. With our simulations we can 
observe all the cases represented in Figures 1 and 2. 

In Figure 3 we show the time evolution of the distribution function p(v, t), in 
the case of a = 1 and b = 0.5 for which there is only one steady state according to 
Theorem 3.1, considering as initial data a Maxwellian with vq — 0 and cTq = 0.25 
in (5.1). We observe that the solution after 3.5 time units numerically achieves the 
steady state with the imposed tolerance. The top left subplot in Figure 4 describes 
the time evolution of the firing rate, which becomes constant after some time. This 
clearly corresponds to the case of a unique locally asymptotically stable stationary 
state. Let us remark that in the right subplot of Figure 3, we can observe the Lipschitz 
behavior of the function at Vr as it should be from the jump in the flux and thus on 
the derivative of the solutions and the stationary states, see Section 3. 

For b — 1.5, we proved in Section 3 that there are two steady states. With our 
simulations we can conjecture that the steady state with larger firing rate is unstable. 
However the stationary solution with low firing rate is locally asymptotically stable. 
We illustrate this situation in the bottom left subplot in Figure 4. Starting with a 




Fig. 3 Distribution functions p(v, t) for b = 0.5 and a = 1 at differents times. 
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Fig. 4 Firing rates N(t) for j = l, Top left: b = 0.5 with initial data a Maxwellian with: vq = 0 and 
Og = 0.25. Top right: b = 3 with initial data a Maxwellian with: ug = — 1 and = 0.5. Bottom left: 
b = 1.5 considering two different initial data: a Maxwellian with: vq = — 1 and er^ = 0.5 and a profile 
given by the expression (3.2) with N = 2.31901. Bottom right: b = —1.5 with initial data a Maxwellian 
with: vq = — 1 and = 0.5. The top right case seems to depict a blow-up phenomena demonstrated in 
Theorem 2.2. 

firing rate close to the high stationary firing value, the solution tends to the low firing 
stationary value. 

In Figure 5 we analyze in more details the behavior of the steady state with larger 
firing rate. The left subplot presents the evolution on time of the firing rate for dif- 
ferent distribution function starting with profiles given by the expression (3.2) with 
N an approximate value of the stationary firing rate. We show that, depending of the 
initial firing rate considered, its behavior is different: tends to the lower steady state 
or goes to infinity. The firing rate for the solution with initial A^o = 2.31901 remains 
almost constant for a period of time. Observe in Figure 5 that the difference between 
the initial data and the distribution function at time f = 1.8 is almost negligible. How- 
ever, the system evolves slowly and at t = 6 the distribution is very close to the the 
lower steady state, see the bottom left subplot in Figure 4. 

In the bottom right subplot of Figure 4 we observe the evolution for a negative 
value of b, where we know that there is always a unique steady state, and its local 
asymptotic stability seems clear from the numerical experiments. 

The number of steady states is related with well-known neuronal phenomena, for 
instance, asynchronous behavior, when there exists only one stationary solution. Let 
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t 



Fig. 5 For b = 1.5 and a{N) = 1 figures show unstability of the steady state with higher firing rate. Left: 
Evolution on time of the firing rate considering different initial firing rate. Right: Evolution on time of the 
distribution function with initial firing rate 2.31901. In both figures we have considered Vr = 1; Vp = 2. 



us mention that bi- and multi-stable networks have been used to describe binocular 
rivalry in visual perception [24] and the process of decision making [25]. Our sim- 
ulations show that the simple NNLIF model (1.4) describes networks with only one 
steady states and several stationary states, in terms of the connectivity parameter b. 

No steady states. - The results in Section 3 indicate that there are no steady states 
for b — 3. In Figure 6 we observe the evolution on time of the distribution function 
p for this choice of the connectivity parameter b. In Figure 4 (right top) we show the 
time evolution of the firing rate, which seems to blow up in finite time. We observe 
how the distribution function becomes more and more picked at Vr and Vp producing 
an increasing value of the firing rate. The synchronization is a phenomenon that is of 
interest to neuroscientists, in these figures we do not observe it, but they do not show 
desynchronization either, since there are not steady states and no convergence of the 
firing rate. Therefore it could be possible to find periodic solutions which describe 
this phenomenon allowing for the formation of point Diracs in the firing rate. In this 
way, it could be related with the situation that we observe in Figures 7 and 8, where 
blow-up is analyzed, since the firing rate looks like tends to be a Dirac Delta. 



0.7 




Fig. 6 Distribution functions p(v, t) for a = 1 and b = 3 at different times. See Figure 4 for the corre- 
sponding plots of N(t). 
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Fig. 7 Parameter values are a = 1 and b = 1.5 and this corresponds to two steady states. Left: Evolution 
of the distribution function p(v, f) in time of an initial Maxwellian centered at no = 1.5 and with variance 
0.005. Right: Time evolution of the firing rate; we observe numerically a blow-up behaviour for an initial 
data satisfying the condition of concentration around Vp described in Theorem 2.2. 

Blow up. - According to our blow-up Theorem 2.2, the blow-up in finite time of 
the solution happens for any value of b > 0 if the initial data is concentrated enough 
on the firing rate. In Figures 7 and 8, we show the evolution on time of the firing rate 
with an initial data with mass concentrated close to Vp for values of b in which there 
are either a unique or two stationary states. The firing rate increases without bound 
up to the computing time. It seems that the blow-up condition in Theorem 2.2 is not 
as restrictive as to say that the initial data is close to a Dirac Delta at Vf, as it is 
showed in Figure 7, where the initial condition is far from Dirac Delta at Vf . Let us 
finally mention that blow-up appears numerically also in case of a (AO = ciq + a\N , 
but here the blow-up scenario is characterized by a break-up of the condition under 
which (1.6) has a unique solution N, that is, 




Fig. 8 Parameter values are a = 1 and b = 0.5 and this corresponds to a single steady state. Left: Evo- 
lution of the distribution function p(v, t) in time of an initial Maxwellian centered at vq = 1.83 and with 
variance 0.003. Right: Time evolution of the firing rate; again this seems to be a typical blow-up behaviour. 
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Therefore, the blow-up in the value of the firing rate appears even if the derivative of 
p at the firing voltage does not diverge. This kind of behavior could be interpreted 
as a synchronization of a part of the network, since the firing rate tends to be a Dirac 
Delta. 



6 Conclusion 

The nonlinear noisy leaky integrate and fire (NNLIF) model is a standard Fokker- 
Planck equation describing spiking events in neuron networks. It was observed nu- 
merically in various places, but never stated as such, that a blow-up phenomena can 
occur in finite time. We have described a class of situations where we can prove 
that this happens. Remarkably, the system can blow-up for all connectivity parameter 
b > 0, whatever is the (stabilizing) noise. 

The nature of this blow-up is not mathematically proved. Nevertheless, our es- 
timates in Lemma 2.3 indicate that it should not come from a vanishing behaviour 
for v s» — co, or a lack of fast decay rate because the second moment in u is con- 
trolled uniformly in blow-up situations. Additionally, numerical evidence is that the 
firing rate N(t) blows-up in finite time whenever a singularity in the system occurs. 
This scenario is compatible with all our theoretical knowledge on the NNLIF and 
in particular with L l estimates on the total network activity (firing rate N(t)). Fur- 
ther understanding of the nature of blow-up behavior, and possible continuation, is a 
challenging mathematical issue. This blow-up phenomenon could be related to syn- 
chronization of the network therefore an interpretation in terms of neurophysiology 
would be interesting. 

On the other hand, we have established that the set of steady states can be empty, 
a single state or two states depending on the network connectivity. These are all com- 
patible with blow-up profile, and when they exist, numerics can exhibit convergence 
to a steady state, that is, to an asynchronous state for the network. Besides better un- 
derstanding of the blow-up phenomena, several questions are left open; is it possible 
to have triple or more steady states? Which of them are stable? Can a bifurcation 
analysis help to understand and compute the set of steady states? 



This appendix is devoted to prove a Hardy-Poincare's-like inequality more gen- 
eral that the one stated in Proposition 4.3. Given m, n > 0 on (0, oo) such that 
f^ c m(y)dy — 1, we want to show that for all functions / on (0, oo), such that 

f^ C m(y)f(y) dy — 0, the following inequality holds: 



provided all integrals make sense. Proposition 4.3 follows by considering m{x) — 
n(x) — K min(x, e~ x Z 2 ), where A" is a constant in such a way that / 0 °° m(y) dy — 1 



Appendix 




(7.1) 
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and parameterizing the interval from (— oo, W] instead of [0, oo). Note that is 
equivalent to the function m both at 0 and at oo in terms of asymptotic behavior. 
We point out that for this kind of functions, m(x) — n{x) — Kmm(x, e~ x I 2 ), the 
Muckenhoupt's criterion for Poincare's inequality or their variants in [32, 33] cannot 
be used, since \/n{x) is not integrable at zero. However, the inequality (7.1) holds 
provided that J 0 °° m(y)f(y)dy — 0. The main ingredient in the proof of (7.1) is to 
write 

2 



= Uf™m(y)(f(x)-f(y))dy\ 



where we used f 0 °°m(y)f{y)dy — 0 and J 0 °° m(y) dy = 1. In this way we can con- 
sider two functions, to be chosen later, cf>, i]/ > 0 to obtain 



^ 2 

dy 



\fW 2 =1(1™ m(y)(f(x) - f(y)) dy 

-{fo m (y) f ff(z) dz dy ) + ( / 00 m(y) L y //(z) dz 

(px pz \ 2 / poo poo 

J f\z)J m(y)dydz) +U f'(z) J m(y)dydz 

< / f(z) z n(z)(p(z)dz / - — - dz 
Jo Jo n(z)(j>(z) 

f°° , , f°° (Pm(y)dy) 2 

+ / f'(z) 2 n(z)if(z)dz / J < r ry 7 v dz. 
Jx Jx n(z)f(z) 



-oo 

m(x)\f(x)\ 2 dx <I + II, 



Therefore, we have 

2 Jo 

where, using the notation, M(x) — m(y) dy, 

\2 



I 

and 



r°° r x . , r* M{ y y 

= / m{x) \ f'{z) 2 n{z)Hz)dz / Ky ' dydx 
Jo Jo Jo n(y)<p(y) 



f 00 f 00 , 9 [°° (1 - M(y)) 2 

11=/ m(x) / f'{z) 2 n{z)ir{z)dz / - ^ sfyJx. 

To conclude the proof we analyse both terms. Considering Fubini's Theorem we ob- 
tain 

1= (°° m{x) f f'(z) 2 n(z)(P(z)dz f* M ^ dydx 
Jo Jo Jo n(y)<p(y) 

f°° , 2 f°° f x M (y) 2 

— / f (z) n(z)(t>{z) / m(x) / ————dydxdz 
Jo Jz Jo n(y)4>(y) 
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with 



and 



with 



rOO 

<Ai / f\z) 2 n(z)dz 
Jo 



f°° M(y) 2 , , 
A\ :=sup0(z) / - — - (1 - M(yvz)Wy where yvz = max(y, z), (7.2) 
z >o Jo "00000 



/•oo poo poo q _ ( v )) 2 

11= / m(x) / f(z) 2 n{z)f{z)dz / ^— -yj^-dydx 
Jo Jx Jx n(y)f(y) 

r t'tV < w n f , , p (1 - M(y)) 2 
= / f (z) n( z )f(z) m(x) — dydxdz 
Jo Jo Jx n(y)f(y) 

POO 

<A 2 / /(z) 2 «(z)Jz 

JO 



f 00 (1 -M(;y)) 2 

A 2 := supi/f(z) / — M(y A z)dy where yAz = min(j, z). (7.3) 

z>o J z n(y)f(y) 

To conclude the proof, it remains to choose the functions <j> and i/r such that, A\, A 2 < 
oo. Using L'Hopital's rule and, after some tedious but easy computations and calculus 

_ 2 

arguments, we obtain that in the case at hand, m(x) = n{x) — K min(x, e~ x ), we can 
take <j>(x) — if(x) — 1 + ■sfx. 
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